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Abstract 

In this paper we describe a version of London Langevin molecular dynamics simulations that 
allows for investigations of the vortex lattice melting transition in the highly anisotropic high- 
temperature superconductor material Bi2Sr2CaCu20g +( 5. We include the full electromagnetic in- 
teraction as well as the Josephson interaction among pancake vortices. We also implement periodic 
boundary conditions in all directions, including the z-axis along which the magnetic field is ap- 
plied. We show how to implement flux cutting and reconnection as an analog to permutations in 
the multilevel Monte Carlo scheme and demonstrate that this process leads to flux entanglement 
that proliferates in the vortex liquid phase. The first-order melting transition of the vortex lattice 
is observed to be in excellent agreement with previous multilevel Monte Carlo simulations. 

PACS numbers: 74.25. Qt, 74.25.Ha, 74.25.Dw, 74.25.Bt 
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I. INTRODUCTION 



There has been a major research effort in recent years to understand the properties of 
high-temperature superconductors which have been discovered during the eighties. High- 
temperature superconductors belong to the class of superconducting materials known as 
type II that allow for partial magnetic flux penetration whenever the external field satisfies 
H c i < H < if ^ 1 ! 2 ! 3 . The flux penetrates the sample in the form of flux-lines (FL's), each 
containing a quantum unit 0o — hc/2e of flux. At low temperature the FL's form an 
ordered hexagonal lattice (Abrikosov lattice) due to their their mutual repulsion. At high 
temperature and/or magnetic field this lattice melts due to thermal fluctuation o 4 i 5 i 6 i 7 i 8 . 

High-temperature superconductors are anisotropic materials which are made from stacks 
of superconducting layers associated with copper-oxide planes. The layers are weekly cou- 
pled to each other. The parameter measuring the anisotropy is 7, defined as 7 2 = m z /m±, 
where m z and m_i_ denote the effective masses of electrons moving along the c axis (perpen- 
dicular to the superconducting planes) and the ah plane, respectively. While for the material 
YBa 2 Cu307_,5 known as YBCO the anisotropy is somewhere between 5-7, for the material 
Bi 2 Sr2CaCu 2 08+5 known as BSCCO, the anisotropy is estimated to be between 10 to a 100 
times larger. 

For BSCCO and highly anisotropic materials similar to it, each FL is represented more 
faithfully by a collection of objects referred to as pancake vortices or just "pancakes"—^. 
Pancakes are centered at the superconducting planes. Each pancake interacts with every 
other pancake, both in the same plane and in different planes. The interaction can be shown 
to consist of two parts. The first part is called the electromagnetic interaction (or simply 
magnetic) and it exists even in the case that the layers of the material are completely 
decoupled, so no current can flow along the c-axis of the sample. The electromagnetic 
interaction originates from screening currents that arise in the same plane were a pancake 
resides as well as in more distant planes. This leads to a repulsive interaction among pancakes 
in the same plane and an attractive interaction among pancakes in different planes*^. 

The second part of the interaction is called the Josephson interaction 2 ! 11 ^ 2 . It results 
from the fact that there is a Josephson current flowing between two superconductors sep- 
arated by an insulator and this current is proportional to the sine of the phase difference 
of the superconducting wave functions. The two superconductors in the present case are 
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the adjacent CuC>2 planes. When two pancakes belonging to the same stack and residing 
in adjacent planes move away from each other, the phase difference that originates causes 
a Josephson current to begin flowing between the planes. This results in an attractive in- 
teraction between pancakes that for distances small compared to r g = jd is approximately 
quadratio^iii in the distance. Here we denoted by d the inter-plane separation and 7 is the 
anisotropy. When the two adjacent pancakes are separated by a distance larger than r g , a 
"Josephson string" is formed, whose energy is proportional to its lengt h 12 ' 13 . 

In three recent paper o 14 i 15 i 16 we presented results of multilevel Monte Carlo simulations 
performed on high temperature superconductors. In the first publication 14 both YBCO 
and BSCCO were treated with and without columnar defects. In that paper we included 
the effect of the electromagnetic interaction only as an in-plane interaction which is valid 
in the approximation that the FL's do not deviate too much from a straight line and the 
anisotropy is not too large. This is usually the case for YBCO 1 ^ and is justified for highly 
anisotropic materials if the anisotropy is not higher than about 250, which is often not the 
case for BSCCO where for optimally doped samples one expects anisotropies in the range 
of 400-500^. 

In the second paper we conducted multilevel Monte Carlo simulations including both 
the in and out of plane electromagnetic interactions plus the Josephson interaction among 
nearest neighbor pancakes in adjacent planes. The Josephson interaction is often neglected 
in simulations of the highly anisotropic BSCCO, but we showed that it is crucial to obtain 
the proper scaling behavior of the results and should not be entirely neglected. We also 
implemented periodic boundary conditions in all directions including the z direction. In 
the third paper a newer approximation to the Josephson coupling has been derived from a 
numerical solution of the two-dimensional sine Gordon equation, which is meant to improve 
on the previous approximation introduced by Ryu et alM. 

Molecular dynamics (MD) is a powerful tool for simulations of physical systems and it 
often serves as an alternative to Monte Carlo (MC) simulations. Its advantage is that it can 
be used to investigate the real dynamics of the system as opposed to MC simulations that 
are used for obtaining equilibrium properties. However MD simulations could be plagued 
by the absence of ergodicity when applied to systems represented by path integrals^ and 
there is also the problem of implementing permutations for the case of identical particles like 
Bosons. The problem of ergodicity is really not much of an issue for Langevin simulations 
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since the thermal noise helps the system explore the configuration space and it can be shown 
by using the corresponding Fokker-Planck equation that equilibrium is reached in the long 
time limit. For flux-lines (FL's) we have found a way to implement "permutations" in the 
MD simulations by flux cutting and recombining as will be explained further below. We 
were also able to implement periodic boundary conditions in all directions (including the z 
direction) and to include the in and out of plane electromagnetic interaction as well as the 
Josephson interaction using the new approximation we have recently obtained 16 . 

Results that clearly show the first-order melting transitions in BSCCO for fields of 100-200 
gauss are presented below. There is an excellent agreement with the results of our multilevel 
Monte Carlo simulations 15 including the proliferation of non-simple loops corresponding to 
flux entanglement above the melting transition. 

At this point we should briefly discuss some previous applications of MD Langevin simu- 
lations for investigations of vortex-lattice phenomena. Wilkins and Jensen 2 ^ used Langevin 
dynamics to investigate the melting transition in the presence of point disorder in layered 
superconductors. However they used the non-realistic Gaussian potential among pancake 
vortices instead of the actual long-ranged logarithmic interactions derived by Clem and 
others*^. They observed a signature of a first order transition that disappear completely 
when the disorder is strong. 

van Otterlo et alM used Langevin dynamics for the case of YBCO where flux-lines rather 
than individual pancakes are the relevant dynamical variables, the electromagnetic interac- 
tion is taken only in plane and then there is the bending forces due to the line tension. The 
authors introduce point disorder and investigate the Bragg glass to vortex glass transition. 

Olson et a/ . 23 ' 24 use Langevin dynamics for pancake vortices in BSCCO. However they 
take into account only the electromagnetic interaction and neglect the Josephson interaction 
entirely, thus effectively using 7 = 00. They also do not implement periodic boundary 
conditions in the z direction, nor do they implement flux-cutting and recombination. Instead 
of varying the magnetic field they use an artificial parameter S m that changes the relative 
strength of the in-plane and out-of-plane interactions. However varying this parameter 
away from unity makes the interaction of a single pancake with a straight stack of pancakes a 
distance R away different from K (R/ X)^. These authors are able to observe the decoupling 
transition of the superconducting planes that occurs at high magnetic fields. They also 
include the effects of point disorder and in addition they investigate the effect of a driving 
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force, like an electric current going through the sample. 

Kolton et alM use MD simulations at T=0 which are therefore not of the Langevin 
type. They implement periodic boundary conditions in all direction and include the full 
long ranged electromagnetic interaction but neglect the Josephson coupling. They study 
current driven pancakes in highly anisotropic superconductors. 

Fangohr et ali^ use both MD Langevin simulations and Monte Carlo to study the melting 
transition in highly anisotropic superconductors. As an alternative to including the full long- 
range electromagnetic interactions they use a mean-field approach in which the instantaneous 
density of pancakes in other layers than the currently simulated layer is replaced by an 
average density, thus leading to an effective "substrate potential"—, that is adjusted self 
consistently. These authors do not include the Josephson interaction. Note that our results 
for the case of infinite anisotropy as discussed in the Results section and in Reflisl agree 
with the results of this paper. 

II. THE MODEL 

The equation of motion for the m'th pancake vortex is 

dr] ^ = _ Vm V({r n }) + fL + Cm{t y pi) 
at 

The pancake label m stands actually for two indices (i,p) where p is the plane label and 
i is the pancake label in that plane. The position R m is a two component vector in the 
plane. Here we have used the over-damped model for vortex motion in which the velocity of 
the vortex is proportional to the applied force and rj is the viscous drag coefficient per unit 
length given by the Bardeen-Stephen^ expression 

>,= ^#, (2.2) 

Pn.C 2 

with p n is the normal-state resistivity, d is the interlayer spacing between CuOi planes that 
is taken to be equal to the width of the pancake vortex. V is the potential energy depending 
on the position of all pancakes and includes both the magnetic energy and Josephson energy. 
The force is minus the gradient of the potential energy with respect to the position of the 
m'th pancake, fi is a driving force (if present), for example the Lorentz force induced by a 
current. C, m is a white thermal noise term which satisfies 

&)&)) = 2kT v d 6 a p6 m J(t - t'). (2.3) 



In Eq. (|2.3|) a and (3 refer to the x and y components of the vector £ and m and n are 
pancake labels, k is Boltzmann's constant. In our simulations we measure distances in units 



of a = y 20o/-By / 3 where I? is the magnetic field. We measure energy in units of end where 
en(T) = (0 o /47tA) 2 is the basic energy scale per unit length and A is the penetration depth. 
We measure time in units of rjal/e . Putting 

R m = a R m ; t = f^A t; V = a^V; V = (e d)V; (2.4) 
Sl = (tod/a )f L ; Cm = ^od/a )C m ; kT = (e d)f. 



we obtain 



with 



<7 Ft 

= -V m V({r n }) + f L + Cm(t). (2.5) 

dt 



<Cm(*)S(?)> = 2T5 aP 5 mn 5(i-t'). (2.6) 



In the simulation we take t to be discreet with an increment At. Thus instead of the Dirac 
delta function 5(t) we take a function which is zero everywhere except when t — 0, in which 
case it is 1/ At. Thus we take 



C(t) = pT/At X a m (t), (2.7) 

where x is a normally distributed random number with zero mean and unit variance. 

To give an example of the magnitude of the various units used we quote their values 
for T = QOK and B = 100G. In that case we have a « 4887 A, e d w 4.685 x l(T 14 erg 
~ 339.5 K/k. Ref. 29] quotes a value for 77 for a single crystal BSCCO of around 1 x 1CT 7 
g/(cm s). Based on this value the time unit is about 0.765ns. The value of the time unit is 
unimportant for the results of the present paper since we report on equilibrium properties. 

We now discuss the expressions used for the various interactions and the methods used 
to implement periodic boundary conditions. 



A. Electromagnetic Coupling 



For the in-plane interaction between two pancakes one has,^^ 



u(fi. J .Q) =21n ^L_^ ln c 

tod Rij A \ Rij 



2 In — - T ( In — - E 1 (R lj ) , (2.8) 
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where Rij = |R» )P — Rj, P | is the radial distance in cylindrical coordinates. Here R is a two 
dimensional vector with components x and y. 

The interaction between two pancakes (Ri )P1 , P\d) and (Rj, P2 , p^d) is given in the case 
when the pancakes are situated at different planes by 



C is some unimportant constant that cancels out upon taking energy differences. We see 
that Ei(Rij) = E2(Rij,0). This form of energy can be derived either by starting from 
Lawrence-Doniach model2*2Jl or by following ClerrA 

We choose our simulation cell to have a rectangular cross section of size ao -J Nji x 



We usually worked with 36 flux lines. The aspect ratio of the cell was chosen to accommo- 
date a triangular lattice without distortion, such that each triangle is equilateral. In the 
z-direction we take N p layers of width d each, where in practice we have chosen N p — 36. 

We now discuss how to implement periodic boundary conditions (PBC) in all directions. 
Let us consider first the implementation of PBC in the z-direction and later we will imple- 
ment PBC in the x and y directions. Periodic boundary conditions mean that every pancake 
interact not only with the actual pancakes in the simulation cell but will all their images in 
other cells which are part of an infinite periodic array. Each image of a pancake is located 
at the same position in the corresponding cell as the original pancake in the simulation cell. 
Thus it is not a reflection through a boundary. 

Let us start with the interaction of a pancake in a certain plane p with another pancake 
in plane p'. Because of the PBC in the z-direction it also interact with all images of p' in 
positions (p' + N p l)d where / is an integer. Thus concerning the first term in Eq. (j2.9J) we 
have to evaluate the sum 




(2.9) 



where = |Rj jPl — Rj, P2 l> an d z = (pi — p-i)d. 

In the above equations we defined the residual interactions 





ao ^J?>Nfi/2 where Nfi is the number of flux lines (number of pancake vortices in each plane). 



oo 



exp(-|Ap|/j) + exp(\Ap\fi) exp(-N p fj,) 
1 — exp(—N p fi) 



f m {Ap)= exp(-|Ap + jV p Z|/z) 



(2.11) 



l=— oo 
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where we put Ap = p — p' and \x = d/X. The dependence of f m on Ap is rather weak. 
For the second term we have to evaluate the sum 

V E 2 (R, (Ap + N p l)d) = / dy/y V exp -^Jy* + (Ap + N p lf (2.12) 

l=-oo jR l d l=-oo V / 

We now make the approximation, valid for the case when dN p <C A that the sum over / can 
be replaced by an integral dl. The estimated error is small for the range of parameters 
under consideration. Recall that for BSCCO, A(T) = Ao/a/(1 — T/T c ) and it is equal to 
3000 — 5000A for the range of temperatures we work with, whereas N p d ~ 500A for the value 
N p = 36 that has been used in the simulations. Changing variables from I to x — Ap + N p l 
we find 

00 i f°° dv r°° 

J2 E 2 (R,(Ap + N p l)d) « — / dxex V (-^x 2 + y 2 ). (2.13) 

Z= _ 0Q 7V P J Rid V J-oo 

We can now change variables from rectangular (x, y) to polar (p, 9) to get 

1 f 77 d6 f°° 
W —J dpexp(-w) = (2.14) 

sin 6) 

2A T /2 dff ( R \ 2A (R\ 

ex P -T— 7? = ^T^o T , (2.15) 



Jq sin 9 \ A sin ^ / c/A'p \ A 



with Kq being the modified Bessel function of second kind of zero's order. The last integral 



jy using the change of variable z — 1/ sin 9 and then referring to formula 
Thus the total contribution to the out-of-plane pair interaction becomes 



was calculated 
3.384/3 in Ref. 



e d A \dN p \XJ J y r ' \R 

It can be checked that to leading order 



2A / _Ae<i 2 



f -w = JvV +ol -Tr ) ) > (2 ' 17) 

since \ Ap\ < N p . As R — > 0, K Q (R/X) f& ln(X/R) + const. Thus requiring that the energy to 
be finite in the limit R — > we replace the prefactor of Kq by f m (Ap), where the difference 
involves only higher order terms, and the correct limits are obtained both when R is small 
and in the limit when R is large and K (R/X) tends to zero. Thus Eq. (j2.16ft is replaced by: 



€$d X V V A y V i? 



We now turn to the interaction of pancakes in the same plane, again concentrating on 
one pancake, and its interaction with another in the same plane and all its images in other 
cells above or below a distance dN p l in the ^-direction. For the images we must use the out 
of plane interaction. Thus apart from the 21n(C/i?) term we have the same calculation as 
above with only difference is that now Ap = 0. Thus we get 

^)~^(c) + * uo) ( Ko («)^(g)). (, 19) 



e d \RJ X J y ' \ "\XJ \R / 

Let us check if we get the right answer for the case of one pancake interacting with a straight 
stack of pancakes a distance R away. Summing all the pair interactions one obtains 

Using the fact that 

£ U&P) « -j, (2.21) 

Ap=0 



U{R) ~ 2K ( *) (2.22) 



one obtains 



e d \ X 

Which is the correct result for the interaction of a pancake with a straight infinite stack see 
Ref. Q Appendix A. 

Consider also a straight stack of pancakes with one pancake from the stack displaced a 
distance R away. The interaction of that pancake with the rest will be in this case 

m^Uf*)-^*)), ( , 23) 



€$d \ \ A J 

to leading order (with correction of order 1/N P ), which coincides with the result obtained 
by Clern^ provided C is chosen appropriately so the energy vanishes as R — ► 0. 

Thus far we only summed over images in the ^-direction. We now have to implement the 
PBC in the transverse direction. In that case 

K (R/X)^G (R/X,L 1 /X), (2.24) 

Where the Green's function Gq(R/X,Lx/X) satisfies London's equation 

(1 - A 2 V 2 )G ( J R/A, VA) = 2nX 2 5{R) (2.25) 



with PBC in the rectangular cell of dimensions L\ x L 2 with L 2 = y / 3-^i/2. Note that Go is 
not spherically symmetric. Similarly in the case of A — > oo one has to replace the logarithm 

by 



HR/C)^G QC {x/L^y/L 2 ) : (2.26) 

which satisfies PBC. Again the expression is derived in the Appendix. In this case an infinite 
constant independent of R has to be subtracted to make the expression finite. Our final 
expression for the pair energy with fully implemented PBC is 

and similarly 

^ - **■ fe< + ( g » (f ■ x) - d'£))- **> 

B. Josephson interaction 

In a recent paper 16 we derived an approximation to the Josephson interaction among 
pancakes in nearest neighbor planes. The approximation is based on a numerical solution of 
the nonlinear sine Gordon equation in two dimensions. A string-like solution corresponding 
to a Josephson string that connects two singularities has been investigated and its energy 
calculated. It is believed that the derived formula constitutes a better approximation to 
the Josephson interaction than the one previously usecU^. The formula obtained for the 
Josephson interaction is 

U Joesephson {R) = e d (1.55 + \n(X/d)) 0.25 (R/r g f In (9r g /R) , R < 2r g 

= e d (1.55 + ln(A/d)) ((R/r g ) - 0.5) , 2r g < R. (2.29) 

where R is the lateral separation of the pancakes and r g = jd where 7 is the anisotropy and 
d is the inter-plane separation. 

Since the Josephson interaction is between nearest neighbor pancakes in adjacent planes 
it is quite straight forward to implement PBC. A pancake at the top plane (N p ) interacts 
with the closest pancake in the bottom plane 1 as well as with a pancake in plane N p — 1. 
When calculating the lateral distance between pancakes we always measure the "shortest 
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distance" denned as follows: If the actual \Ax\ separation is larger than Li/2 we subtract 
or add L\ depending on the sign of Ax, and similarly for Ay with L2 replacing L\. This 
way the correct distance is obtained even when the adjacent pancake in the plane above has 
exited the simulation cell and emerged close to the other side of the simulation cell. This 
is because when a pancake exits the cell from one side it (or what was its image) enters 
the cell from the other side, so for the Josephson interaction pancakes close to two distance 
boundaries can actually be neighbors. 

In the case of R 3> r g , string-string interactions that involve three and four-body inter- 
actions become important^. However near the melting transition for the range of magnetic 
fields investigated in this paper R « 0.25a ~ 1000A whereas r g 5625A. Thus large trans- 
verse fluctuations for which the string-string interactions become important are statistically 
rare and can be neglected. 

III. DETAILS OF THE SIMULATIONS 

The simulation cell is divided into a 800 x 692 mesh of small cells of area hxh each where 
h = ^Nji/800 (in units of a ), and we tabulate the functions Go and Gqc m each small cell 
thus creating two large 800 x 692 matrices. During the simulations we use the tables as a 
lookup to calculate the pair interaction. For each table we also calculate the negative of the 
gradient and save the two components of the gradient in their own tables. We also tabulate 
the Josephson interaction and its gradient. When simulating we allow pancakes to move to 
arbitrary real locations but in order to calculate the forces we divide the actual position by 
h and round to the nearest integer to use for the lookup tables. 

In each simulation step we move all the pancakes at the same time, using the instan- 
taneous forces. This is done using a time step At. It is very important to chose the time 
step correctly. Consider the magnitude of the white thermal noise given in Eq. (j2.7j) . On the 
average, the distance a pancake moves during a time At is given by v / 2TAt. We choose 
this distance to be either 5h or 8h as explained below. We have used two methods: First we 
have employed a simple Euler method. For a given configuration of pancake we calculate 
the force on each pancake due to the pair potential due all other pancakes, both magnetic 
and Josephson. To this force we add the constant driving force (if any) and the thermal 
noise. Based on this forces we move each and every pancake simultaneously (in parallel) by 
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a distance given by the total force acting on it just before the move times At = (5h) 2 / (2T). 
Then we calculate the new forces and repeat. At each step we generate the thermal noise by 
calling a gaussian vector random number generator to obtain a vector of length 2 x Nfi x N p 
filled with random numbers (total number of pancakes times two force components). The 
Euler method works adequately but is relatively slow since in order to get good results we 
needed to simulate up to a time span of 72 time units or more. 

We found that we can improve performance by using a second order Runge-Kutta method. 
We have tried to use a forth-order method but there has been no further improvement over 
the second order method for reasons that will be explained below. In the second order 
method we use a larger time step of length At = {8h) 2 / (2T), which is about 2.6 larger than 
in the Euler method. In this method we first consider a virtual move of duration At/2 with 
the initial instantaneous forces, we then calculate the new forces at the end of the virtual 
move and we use these forces to move a full time step starting at the original initial position. 
The random noise is only generated once at the original point. Exactly the same random 
noise is used both for the virtual move and for the actual move. Thus the vector of random 
numbers is saved and used in the same order for the virtual and actual moves. Of course 
each move now takes about twice the cpu time than before, but we gain both because the 
results are more accurate and because we use a larger time step that reduces cpu time for 
the same total time span. In this method we could get reliable results in about half to two 
thirds the cpu time needed for Euler's method. 

We believe that the reason we did not get an improvement with the fourth order method 
is that in order to get a reduction in cpu time we need to increase the time step by at least 
a factor of two since each update move requires four evaluations of the forces instead of two. 
But for the same total time span this reduces the number of steps by at least a factor of 
two. This interferes with the statistics of averaging over the thermal noise since the number 
of steps is not large enough to get good statistics so the results are actually not as good as 
the results from the second order method (for the same total time span). 

It is a nice feature of multilevel Monte Carlo that one can implement flux cutting and 
permutations^*^^^, so that flux lines with PBC in the z direction don't end on themselves 
but form loops that wind more than once across the system. We term such loops non- 
simple or "composite loops". For Bosons, the abundance of such loops characterizes the 
superfluid phase 20 . They represent permutations of the particles that differ from the identity 
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permutation. There is an approximate mapping from the world lines of bosons propagating 
along the Euclidean time direction to FL's stretching along the z directio n 17 ! 31 . For flux-lines 
these composite loops represent the entangled state of the vortex liquid above the melting 
transition. In order for this concept to exist one must not neglect the Josephson interactions 
even for a highly anisotropic material like BSCCO since it is the Josephson interactions that 
really tie up a stack of independent pancakes into a flux-line even if loosely so since this 
interaction is weak. 

It is sometimes argued that one can not implement permutations in a molecular dynamic 
simulation as the motion through permutation space is discrete and molecular dynamics 
involves continuous evolution. However in our situation it is quite possible to introduce 
"permutations" in our system and see that they proliferate above the melting transitions. 
In fact the results obtained for the number of loops not ending on themselves agree amazingly 
well with the corresponding results from our Monte Carlo simulations of the same system. 

The way we implement "permutations" is through flux cutting and recombination. We 
assume that within the coupled-planes model, vortices may switch connections to lower 
their elastic energy (in this case Josephson energy) when they cross each other^-. In the 
simulation we construct two matrices of size Nfi x N p which we call the "up" matrix and 
the "down" matrix. For a given pancake % in plane p the "up" matrix points to the pancake 
in the plane p + 1 (or 1 if p — N p ) that is connected to the given pancake (i,p) via a 
Josephson interaction. Generally this is the pancake closest to the given pancake in the next 
plane. The "down" matrix similarly points to the closest pancake below (or in plane N p for 
p — 1). When we start from an initial configuration in which the FL are a straight stack 
of pancake the matrices simply point to the pancake just above or below a given pancake. 
When constructing the force matrix after each time step we check if indeed the "up" matrix 
points to the closest pancake above. If there is a closer pancake than the one given by the 
pointer then we find out its parent in the plane p by using the down matrix, and we check if 
switching the two connections will decrease the sum of the squares of the two distances. If 
it does we cut and switch connections and update the "up" and "down" matrices. We term 
this precess an "exchange" . The reason we use the square of the distances is that in most 
instances the Josephson interaction is proportional to the square of the transverse distance 
(see Eq. (j2.29J) above). This procedure mimics the actual dynamics in which we expect the 
magnetic flux to choose a path that minimizes the Josephson energy. We implement the 
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flux cutting procedure after every update move of the system, but not during the virtual 
half-step in the Runge-Kutta procedure. 

Note that the extent that flux cutting and reconnecting occurs in real experimental 
samples and the existence of the entangled state is still a debatable issued. Flux cutting can 
allow an entangled state to disentangle and vice versa. Our simulations show that just below 
the melting transition, even though some exchanges occur, they soon reverse themselves in 
space or in time, and thus they do not lead to what we refer to as an entangled state 
where composite loops or permutations are abundant. On the other hand when exchanges 
proliferate through the system, a phenomenon that occurs in our simulations just above the 
melting transition, the exchanges do not undo each other, and the system of FL's changes 
from being composed of simple loops each made up of a single FL, to a system composed 
mainly of composite loops that wind up several times around the simulation cell in the z 
direction before returning to the original point. The reason that the exchanges proliferate 
above the melting transition is that the transverse fluctuations become strong enough to 
overcome the potential barriers due the repulsion among pancakes residing in the same 
plane and thus the crossing of FL's occur. 

Crabtree and Nelson^ give a rough back of the envelope estimate of the magnetic fields 
B x \ and B X 2 such that when B x \ < B < B X 2 entanglement should occur in the flux liquid 
phase. For the system size and the values of parameters and field range that we use (100G — 
300G), we verified that indeed the FL liquid should indeed be entangled. It should be 
noted that Wilkin and Jensen^i measured flux cutting by simply observing the rate that 
the nearest neighbors of a given pancake in adjacent planes change during the course of a 
given time interval. In the liquid state many of those events occurred. However they did 
not implement "exchanges", nor did they keep track of composite loops and their relative 
abundance compared to simple loops as we do in our simulations. 

IV. MEASURED QUANTITIES 

We measured the following physical quantities. For details the reader is referred to our 
earlier work 14,15 . 
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A. Energy 



The average energy was obtained by adding the electromagnetic energy of all pairs of 
pancakes combined with the Josephson energy of nearest neighbor pancakes in adjacent 
planes. 

B. Translational structure factor 

The translational structure factors S'(Qj) is defined as, 

S{Qi) = _L_ e (*Q,(R^-R fc , P )) \ (41) 

p f l \jk,p I 
where (...) stands for the time average, and Qj, i — 1,2 stand for the basic reciprocal 
lattices vectors which are given by 

= r-TE ( e i ~ e j cos 8) » ( 4 - 2 ) 

a sin v 

where i,j = (1,2) or (2, 1), 6 = tt/3, a is the size of the unit cell of the triangular lattice 
and ei j2 are the unit vectors along the rhombic unit cell such that 

ei • e 2 = cos#. (4.3) 

Notice that we normalized the structure factor to unity instead of Nfi. We actually try 
different orientations of e.\ to allow for situations that the lattice unit cell does not align with 
the simulation cell and numerically find the angle for which the average (S(Qi) + S{Q2))/2 
is maximal. We then record this value as the measure of translational order. 

C. Mean square deviations 

For each individual flux-line we define the position of the lateral center of mass as Rcm = 
R(i,p)/N p where the sum goes over all the pancake belonging to it. We then define the 
mean square deviations as 

R) = (J2(R^p) ~ Rcm) 2 ) , (4.4) 
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Where the sum is over all pancakes belonging to an individual flux-line and the average is 
over all flux-lines of the system and then taking a time average. The melting transition is 
expected to occur when this quantity satisfies 

Rf/ao > c L , (4.5) 

where Cl is the Lindemann coefficient. 

D. Line entanglement 

As we allow for flux cutting and recombination, we can define the number N e /Nfi as that 
fraction of the total number of FL's which belong to loops that are bigger than the size of a 
"simple" loop. A simple loop is defined as a set of N p beads connected end to end (due to 
the periodic boundary conditions in the z direction), N p being the total number of planes. 
Loops of size 2N p , 3N p ... start proliferating at and above the melting temperature. 

E. Parameters 

Parameters for BSCCO were taken as follows: Ao = 1700 A, d — 15 A and T c = 90 K. 
The temperature dependence of A in this work was taken to follow the Ginzburg-Landau 
convention A 2 (T) = Aq/(1 — T/T c ). See discussion in Ref. 115] 

on the agreement of this 
choice with experiments. For the anisotropy we have used values of 250-400. 

V. RESULTS 

In this section we display some of the results for the melting transition obtained with 
the molecular dynamics method. The case of B — 100 gauss and 7 = 400 is depicted in 
Fig. 1. In subfigure (a) we see the decay of the normalized structure factor. The melting 
temperature is about 69K (corresponding to a reduced temperature of 300K). In subfigure 
(b) we observe the quantity Rj defined above that measures the square of the transverse 
deviations from a straight line. We see that at the transition the Lindemann parameter cl 
is about 0.25 (its square is about 0.06). In subfigure (c) we observe that composite loops 
corresponding to line entanglement start to proliferate above the melting transitions. In 
part (d) we show the jump in the Josephson energy corresponding to a first order transition. 
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FIG. 1: Results for 7 = 400 and B = 100 G. The following quantities are shown: (a) the trans- 
lational structure factor normalized to unity, (b) the mean square transverse deviations about the 
FL's center of mass in units of Oq, (c) the fraction of composite loops as a measure of FL entangle- 
ment, (d) Josephson energy per pancake in units of eod. The temperature is measured in Kelvin 
and so is the reduced temperature T/(l — T/T c ). 

There is a corresponding jump in the total energy that is more difficult to observe since it 
is fractionally smaller. The jumps are of course smoothened by fine size effects, i.e. the fact 
that we have 36 FL's and 36 planes for a total of 1296 pancake vortices. 

The results agree with multilevel MC simulations carried by us and the fact that the frac- 
tion of non-simple loops agrees with the MC shows that "permutations" were implementec 
faithfully in the MD simulations. Some of our previous MC results are given in Ref. ^5. 
Note that in that paper we used a different approximation for the Josephson interaction as 
given by Ryu at al. and hence one has to adjust the values of the anisotropies in that paper 
by about 1.5 to correspond to the current simulation which treats the Josephson interaction 
according to the approximation given in Ref. Q. More recent MC are presented in Ref. 34 
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FIG. 2: (color online) Hysteresis loop displaying the structure factor for 7 = 400 and B = 100 G. 
The direction of the heating and cooling cycle is indicated by arrows. The temperature is changed 
by 0.5K increments. 

which uses the current scheme. 

In the MC simulations we investigated the finite size effects in more detail. We tried 
to increase the number of FLs to 64 instead of 36. We also simulated with the number of 
planes equal 25, 36 and 50. As the number of FLs and planes increase the transition becomes 
sharper but its position does not move by more than IK from its value for 36 FLs and 36 
planes which we use in the current simulations. Our aim here is not to pinpoint the melting 
transition to a high accuracy but to have a simulation method that gives reasonable results, 
and can be the basis for simulations on larger systems if one needs to obtain better precision. 
The equilibration times in the MD simulation were chosen to give a good agreement with 
the MC simulations. We also observed that if the equilibration time is not long enough the 
melting appears gradual. By increasing the simulation time the transition becomes sharper 
up to a point when increasing the equilibration time further has no noticeable effect on the 
results. That is how we fixed the equilibration time. Usually it corresponds to at least 
10,000 MD moves (in each move all the pancakes are moved at once) out of which 5000 
moves are discarded before the measurement process begins. 

Since the melting transition is a first order transition we expect hysteresis effects if we 
perform a heating an cooling cycle. The hysteresis should be enhanced by the fact that the 
FLs in the liquid phase are entangled and it takes considerable time for them to disentangle. 
Most of our simulations were done on a parallel machine where at each temperature we start 
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FIG. 3: Results for 7 = 400 and B = 200 G. The same quantities are shown as in Fig.l. 

from an ordered vortex configuration. However, in order to observe the hysteresis we carried 
out a heating and cooling cycle at 0.5K increments where at each temperature we started 
from the last configuration obtained in the previous temperature. We simulated for 72 time 
units at each temperature. The results for B = 100 G and 7 = 400 are depicted in Fig. 2. 

In Fig. 3 we see the melting transition for B = 200G and as expected it occurs at lower 
temperature. From the figure one can read an approximate transition temperature of 
(corresponding to a reduced temperature of 210-fT). To simulate each point in the above 
figures took between 12-24 processor-hours on a 1GHZ processor. Time spans were between 
72-108 time units (in the units discussed in Sec. 2 above), and about half of this time was 
discarded for equilibration and half used for measurements. The cpu times are larger by 
of factor of 2-3 compared with the corresponding times in our multilevel MC simulations 
because of the need to calculate all of the forces, not just the energies, however since 
this method can be used to implement real dynamics in addition to the measurement of 
equilibrium properties only as done in MC simulations, the extra time can certainly be 



19 




FIG. 4: (color online) Phase diagram showing the MD simulation results for 7 = 400 (circles), the 
MD simulation results for 7 = 00 with no Josephson coupling (triangles) and the experimental 



results (squares) of Ref. 
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tolerated. 

In Fig. 4 we display the simulation results for 7 = 400 and 7 = 00 (no Josephson cou- 
pling) as compared to the experimental results of Majer et a/.—. For a more complete phase 
diagram obtained using the MC method for different values of the anisotropy parameter 
and more values of the -B-field see Ref. [lJ5. Notice that although we have chosen 7 = 400 
in order that the melting temperature for B = 150 G will roughly agree with experimental 
result o^ 1 ^ , the simulated melting curve is steeper than the observed experimental melting 
curve in pristine systems. This has been observed and discussed before^. We should re- 
member that experimental pristine systems always include a certain amount of point defects 
that tend to reduce the melting temperature. The effectiveness of these defects increases 
when the temperature is decreased and this causes the melting curve to flatten down in the 
experimental curves of the phase boundary in the B-T plane as compared with the the- 
oretical results for a defect free system. The experimental "irreversibility line" which lies 
just below the melting line is steeper and agrees better with the simulations. In Ref. Q 
we also showed that when the Josephson interaction is present the data for the melting line 
for different anisotropies collapses onto a single straight line when ln(i?7 2 ) is plotted versus 
ln(kT / 6qcI) . This confirms a prediction of Koshelev^I that when the Josephson interaction is 
important the phase boundary is given by a single dimensionless function of the dimension- 
less parameters {kT/e^d) and r g /al oc B^ 2 . For 7 = 00 our results are in agreement with 
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Dodgson et alM and scaling is obtained when plotting BX 2 /cf) vs kT/e d. 

VI. CONCLUSIONS 

In this paper we have shown that molecular dynamics is a powerful tool that can be used to 
obtain the properties of the melting transition in a similar way to multilevel MC simulations. 
We showed how to implement flux cutting and recombination and obtained results showing 
flux-line entanglement similar to those obtained by implementing permutations in the MC 
simulations. We have included both the electromagnetic interaction among all pancakes and 
the Josephson interaction among nearest neighbor pancakes in adjacent planes. We have 
implemented periodic boundary conditions in all directions. 

Our next goal is to include defects, either in the form of columnar defects and/or point 
defects and to investigate steady-state, non-equilibrium properties of the system when a 
current is flowing. 
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APPENDIX A: ENERGY SUM OVER THE IMAGES 



Here unlike our previous papers we work with a rectangular simulation cell with edges of 
size Li and L 2 = Li\f?>/2. The function Go is a solution to the London equation 

(1 - A 2 V 2 )G (R, A) = 2vrA 2 5(R), (Al) 

with the parameter A (penetration depth) setting the scale for the range of the interaction. 
Periodic boundary conditions are to be satisfied in the x and y directions. The solution is 
given by 

2ttA 2 exp(«Q ■ R) 
Q 

where Q = ni(2Tr/Li)i + n 2 (27i/L 2 )j is a reciprocal lattice vector and i%\ and n 2 are integers. 
The summation over rii can be done analytically using a well known formula (see Gradshteyn 
and Ryzhik^, Eq.(l. 445/2)). We are left with one summation: 

n it) w L i ^ cosh(a n (7r - t x )) cos(nt 2 ) L x cosh(a (vr - t x )) 

G (K, A) = — > — r V— — 7—, (Ad) 

L 2 ^— ^ a n smh(a n 7r) 2L 2 a smn(a 7r) 

where we defined 



W)^E^^, (A2) 



L\ / „ Lo 2nx 2ny , . . 

a " = ^V" 2 + i^' '■ = — • i2 = T? (A4) 

and < x < Li, < y < L 2 . We used this formula, and a similar one obtained by first 
summing over n 2 , to calculate Go numerically for finite A. In the limit A — > oo we can obtain 
an equation for the "periodic logarithm". In that limit we have a n = nL 1 /L 2 but we see 
that the last term in Eq. (jA3|) diverges. The diverging term L\j (L 2 2nal) is independent of 
the position R and can be subtracted out. The final expression for Goc* is 

x y Li ^ cosh(a n (7r - ti)) cos(nt 2 ) , vrLi / 3ti 3tj\ 
Li L 2 L 2 ' a n smh(a n 7r) 6L 2 \ n 2% 2 I 

n=l • ' 



with a n = nLi/L 2 . 
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